# ******** Regressions in "Substitution Between Money and Near-Money Assets". *******
# Created by Wenhao Li, 2021 CC.
# Updated on March 26, 2022.  Ready for RFS Revision Round 2.
# =================== Load Packages ===================
remove(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("data_definition_header.R")

# Limit the data period to the end of 2019.
TB = TB[Year<=2019 & !is.na(Deposits_Total)] 


# ********** Plots ***********
# ===================== Figures: Time Series Plots for Main Variables =====================
# This figure shows the variations in the explanatory variables in the paper

# ===================== Figure 1: Treasury Liquidity Premium =====================
pdf( paste0(FigurePath, "liq_prem_main.pdf"), width=7.5, height =4.5 )
par(mar = c(2, 3.5, 1, 1), mfrow=c(1,1))
MyPlot( TB_full$YearMon, TB_full[, list(LiqPrem) ],  ylab = "Liquidity Premium (%)" )
abline( v = as.Date(as.yearmon( c(WWII_periods[1], WWII_periods[length(WWII_periods)])  )), lty=2  )
text( as.Date(as.yearmon(mean(WWII_periods))), 2, "WWII" )
dev.off()

# ===================== Figure 2: Treasury Debt and Financial Sector Deposits =====================
pdf( paste0(FigurePath, "deposit_debt_quantities.pdf"), width=7.5, height = 4.5 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB$YearMon, TB[, list(Deposits_GDP, DebtToGDP_raw  , DebtToGDP, KVJ_fin_sector_debt_to_GDP)], LegendNames = c( "Deposits/GDP", "Total Treasury Debt/GDP","Net Treasury Debt/GDP", "KVJ Deposits/GDP" ), 
        abline_v = as.Date(as.yearmon( c(WWII_periods[1], WWII_periods[length(WWII_periods)])  )), LegendPosition = "topright" )
text( as.Date(as.yearmon(mean(WWII_periods))), 0.2, "WWII" )
dev.off()

pdf( paste0(FigurePath, "three_quantities.pdf"), width=7.5, height = 4.5 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB$YearMon, TB[, list(Deposits_GDP, DebtToGDP_raw, pmax(KVJ_net_deposits,0))], LegendNames = c( "Deposits/GDP", "Total Treasury Debt/GDP","Non-Bank Liquid Liabilities/GDP"), 
        abline_v = as.Date(as.yearmon( c(WWII_periods[1], WWII_periods[length(WWII_periods)])  )), LegendPosition = "topright" )
text( as.Date(as.yearmon(mean(WWII_periods))), 0.2, "WWII" )
dev.off()


# ===================== Figure 3: Comparison of Deposits Rates and Volume =====================

pdf( paste0(FigurePath, "DepositsComposition.pdf"), width=4.5, height =4 )
par(mar = c(2, 2, 0.5, 0.5), mfrow=c(1,1))
MyPlot( TB$YearMon,   TB[ , list( checking/GDP_nominal, saving/GDP_nominal, time/GDP_nominal )  ], LegendNames = c( "checking deposits/GDP", "savings deposits/GDP", "small time deposits/GDP" ) )
abline( v = as.Date(as.yearmon( c(WWII_periods[1], WWII_periods[length(WWII_periods)])  )), lty=2  )
text( as.Date(as.yearmon(mean(WWII_periods))), 0.4, "WWII" )
dev.off()
pdf( paste0(FigurePath, "DepositsRates.pdf"), width=4.5, height =4 )
par(mar = c(2, 2, 0.5, 0.5), mfrow=c(1,1))
MyPlot( TB$YearMon, TB[ , list( checking_rate, saving_rate, smalltime_rate , FFR) ], LegendNames = c( "checking deposit rate", "savings deposit rate", "small time deposit rate" , "Fed Funds rate" ), colors=c("blue", "red", "purple", "black"), ylim=c(-0.5, 13) )
dev.off()

#   Not appear in the paper
pdf( paste0(FigurePath, "Tbill_vs_Devt.pdf"), width=6, height =4 )
par(mar = c(2, 3.5, 0.5, 0.5), mfrow=c(1,1))
MyPlot( TB$YearMon, TB[ , list( DebtToDeposits_standardized, scale(TbillToDeposits)) ], LegendNames = c( "Debt/Deposits", "Tbill/Deposits" ),ylab= "Standardized Deviations")
dev.off()


# ===================== TABLE 1: Log Regressions =====================
index =  !is.element(TB$Year,  WWII_periods) | is.element(TB$Year,  WWII_periods) 
index = index & is.element(1:nrow(TB), seq(3,nrow(TB),by=3))   # The second one picks last month in each quarter
# Note: index now is the whole sample, but at quarterly frequency.
regs=list()
regs[[1]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX)  ,  data=TB[index]  )
regs[[2]] = lm(  LiqPrem_Log  ~  log(VIX) + log(DebtToDeposits) ,  data=TB[index]  )
regs[[3]] = lm(  LiqPrem_Log  ~  log(VIX) + log(Debt_to_KVJ_Deposits) ,  data=TB[index]  )
regs[[4]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX) + log(DebtToDeposits) ,  data=TB[index]  )
regs[[5]] = lm(  LiqPrem_Log  ~  log(FFR) + log(VIX) + log(Debt_to_KVJ_Deposits) ,  data=TB[index]  )

vars.order = c("log\\(DebtToDeposits\\)","log\\(Debt_to_KVJ_Deposits\\)","log\\(FFR\\)","log\\(VIX\\)","Constant")
stargazer( regs, omit = c("f"), omit.stat = c("F", "adj.rsq", "ser"), omit.table.layout="n",
           se=RobustSE(regs), df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0), 
           covariate.labels = c(  "log($ \\frac{\\text{Net Tsy}_t}{\\text{Deposits}_t} $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{KVJ Deposits}_t} $)","log(FFR$_t$)", "log(VIX$_t $)" ), 
           dep.var.caption = "\\textit{Dependent variable:} log(Treasury liquidity premium$_t$)", dep.var.labels.include = FALSE , 
           order = paste0("^",vars.order,"$") )


# ===================== TABLE 2: OLS and IV Regressions in Logs =====================
regs_OLS=list()
regs_OLS[[1]] = lm(  LiqPrem_Log  ~  log(FFR)  + log(DebtToDeposits) + log(VIX) ,  data=TB[index]  )
regs_OLS[[2]] = lm(  LiqPrem_Log  ~  log(FFR)  + log(Debt_to_KVJ_Deposits) + log(VIX) ,  data=TB[index]  )

IV_formula_level = "log(FFR) + log(DebtToGDP_raw) + DebtToGDP_raw  + log(VIX)"
regs = list()
regs_monthly_dummies = list()
IV_formulas = list()
quantity_ratios = c("log(DebtToDeposits)",  "log(Debt_to_KVJ_Deposits)" )
for(i in 1:length(quantity_ratios)){
  IV_formulas[[i]] = paste0( "LiqPrem_Log ~ log(FFR) + ", quantity_ratios[[i]], " + log(VIX)", " | ",  IV_formula_level)
  regs[[i]] = ivreg( IV_formulas[[i]] , data = TB[index])
}
# Output in text
vars.order = c("log\\(DebtToDeposits\\)","log\\(Debt_to_KVJ_Deposits\\)","log\\(FFR\\)","log\\(VIX\\)","Constant")

stargazer( c(regs_OLS, regs), type="text", se=c(RobustSE(regs_OLS), RobustSE(regs)), 
           df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0) , omit.stat = c("F", "adj.rsq", "ser"), 
           dep.var.caption = "\\textit{Dependent variable:} log(liquidity premium$_t$)", dep.var.labels.include = FALSE, no.space = TRUE,
           order = paste0("^",vars.order,"$") )

stargazer( c(regs_OLS, regs), se=c(RobustSE(regs_OLS), RobustSE(regs)), 
           #  c( "~~Cragg-Donald statistic", c( rep("",length(regs_OLS)) , CD_stats_vec_from_ivregs(regs) ) ), 
           df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0) , omit.stat = c("F", "adj.rsq", "ser"), omit.table.layout="n", 
           covariate.labels = c(  "log($ \\frac{\\text{Net Tsy}_t}{\\text{Deposits}_t} $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{KVJ Deposits}_t} $)","log(FFR$_t$)", "log(VIX$_t $)" ), 
           dep.var.caption = "\\textit{Dependent variable:} log(liquidity premium$_t$)", dep.var.labels.include = FALSE,
           order = paste0("^",vars.order,"$"))


## ************************* Instrumental Variable Regressions and Difference Regressions *************************
## This part implements the same instrumental regressions as Nagel for comparison. 
# **** Fed Funds Futures Data
TB_FF_Futures = as.data.table(read_xlsx("Fed_Fund_Futures(Bloomberg).xlsx" ))
TB_FF_Futures[ , FF1:=100-`FF1 Comdty` ]
TB_FF_Futures[ , FF2:=100-`FF2 Comdty` ]
TB_FF_Futures[ , FF3:=100-`FF3 Comdty` ]
TB_FF_Futures[ , FF4:=100-`FF4 Comdty` ]
TB_FF_Futures[ , YearMon:=as.yearmon(Date) ]
TB_FF_Futures_monthly = TB_FF_Futures[ ,  lapply(.SD,  mean_fun  )   ,by=YearMon ]    # Note: it is critical how we aggregate the Fed funds futures data into monthly data. 
TB_FF_Futures_monthly[ , Date:=NULL ]
TB_FF_Futures_monthly[ ,YearMon := as.Date(YearMon) ]
TB_new = merge( TB,  TB_FF_Futures_monthly, by="YearMon",  all.x=TRUE )

# *** Define the difference variables
TB_Diff = TB_new[ , lapply(.SD, mean_fun) , by=YearMon ]
TB_Diff[ , LiqPrem_Diff:= 100*GenDiff(LiqPrem)   ]
TB_Diff[ , LiqPrem_LogDiff:= GenDiff(LiqPrem_Log)   ]  
TB_Diff[ , LiqPrem_LogWinsorized_Diff:= GenLogDiff(LiqPrem_winsorized)   ]  
TB_Diff[ , DebtToGDP_LogDiff:= GenLogDiff(DebtToGDP)   ]
TB_Diff[ , DebtToDeposits_LogDiff:= GenLogDiff(DebtToDeposits)   ]
TB_Diff[ , DebtToDeposits_LogDiff_Lag:= GenLogDiff(DebtToDeposits,-1)   ]
TB_Diff[ , FFR_Diff:= GenDiff(FFR)   ]
TB_Diff[ , FFR_LogDiff:= GenLogDiff(FFR)   ]
TB_Diff[ , VIX_Diff:= GenDiff(VIX)   ]
TB_Diff[ , VIX_LogDiff:= GenLogDiff(VIX)   ]
TB_Diff[ , TbillToGDP_LogDiff:= GenLogDiff(TbillToGDP)   ]
TB_Diff[ , TbillToDeposits_LogDiff:= GenLogDiff(TbillToDeposits)   ]
TB_Diff[ , TbillToGDP_LogDiff_Lag:= GenLogDiff(TbillToGDP,-1)   ]  # values a month ago
TB_Diff[ , TbillToDeposits_LogDiff_Lag:= GenLogDiff(TbillToDeposits,-1)   ]  # values a month ago

TB_Diff[ , TbillTo_KVJDeposits:= TbillToGDP/KVJ_fin_sector_debt_to_GDP   ]
TB_Diff[ , TbillTo_KVJDeposits_LogDiff:= GenLogDiff(TbillTo_KVJDeposits)   ]
TB_Diff[ , DebtTo_KVJDeposits:= Debt_to_KVJ_Deposits   ]
TB_Diff[ , DebtTo_KVJDeposits_LogDiff:= GenLogDiff(DebtTo_KVJDeposits)   ]
TB_Diff[ , DebtToGDP_raw_LogDiff := GenLogDiff(DebtToGDP_raw) ]

TB_Diff[ , YearMon:=as.Date(YearMon)]
TB_Diff[ , DummyIndex:=month(as.POSIXlt(YearMon, format="%d/%m/%Y")) ]

TB_Diff[ , FF_implied_change := FF3 - FF2 ]
TB_Diff[ , FF_implied_change := c( NA, FF_implied_change[1:(nrow(TB_Diff)-1)] ) ]

TB_Diff[ , FF_Futures_IV_Level :=  c( NA, FF2[1:(nrow(TB_Diff)-1)] ) ]
# TB_Diff[ , FFfutures_IV := c(NA, FF_implied_change[1:(nrow(TB_Diff)-1)]) ]
TB_Diff[ , FFfutures_IV := FF_implied_change ]
TB_Diff[ , FF_Futures_IV_Level_lag_3M :=  c( rep(NA,3), FF4[1:(nrow(TB_Diff)-3)] ) ]


# *** Regression settings with IVs (monthly dummies)
name_list = paste0( "M", seq(1,11) )
IV_formula = name_list[1]
for( i in 1:length(name_list) ){
  TB_Diff[ , (name_list[i]):= DummyIndex==i]
  if( i>=2 ){
    IV_formula = paste0(IV_formula, " + ", name_list[i])
  }
}
IV_formula_monthly_dummies = IV_formula
IV_formula = paste0( IV_formula, " + FFfutures_IV" )

# ===================== TABLE 3: Difference Regressions: OLS and IV ======
# *** Limit to subsamples for comparison with Nagel(2016) on first stage
TB_repoT = TB_Diff[ YearMon >= as.Date("1991-05-01") & YearMon <= as.Date("2020-11-01")  ]  # Note: make the dates fully accurate
index = TB_repoT$YearMon >= as.Date("1991-06-01")   # Note: replicate Nagel's date range. 

regs = list()
regs[[1]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff , data = TB_repoT)
regs[[2]] = lm( LiqPrem_Diff ~  VIX_Diff +  TbillToDeposits_LogDiff , data = TB_repoT)
regs[[3]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff +  TbillToDeposits_LogDiff , data = TB_repoT)
regs[[4]] = lm( LiqPrem_Diff ~  FFR_Diff + VIX_Diff +  TbillToDeposits_LogDiff + TbillToDeposits_LogDiff_Lag , data = TB_repoT[index])
IV_formula1 = paste0( "LiqPrem_Diff ~ FFR_Diff | ",  IV_formula)
IV_formula2 = paste0( "LiqPrem_Diff ~ TbillToDeposits_LogDiff | ",  IV_formula)
IV_formula3 = paste0( "LiqPrem_Diff ~ FFR_Diff + TbillToDeposits_LogDiff | ",  IV_formula)
IV_formula4 = paste0( "LiqPrem_Diff ~ FFR_Diff + TbillToDeposits_LogDiff + TbillToDeposits_LogDiff_Lag | ",  IV_formula)
regs[[5]] = ivreg( IV_formula1 , data = TB_repoT[index])
regs[[6]] = ivreg( IV_formula2 , data = TB_repoT[index])
regs[[7]] = ivreg( IV_formula3 , data = TB_repoT[index])
regs[[8]] = ivreg( IV_formula4 , data = TB_repoT[index])

vars.order = c("TbillToDeposits_LogDiff","TbillToDeposits_LogDiff_Lag","FFR_Diff","VIX_Diff")
stargazer( regs,  se=RobustSE(regs,lag_input = 12) , add.lines = list(  c("Weak instruments test", rep("",8)),  c( "~~CD statistic", CD_stats_vec_from_ivregs(regs) ),  c( "~~Critical value", rep("",4), 6.56 , 6.56, 6.23, 5.87 )  ),   # c("Instruments", rep("",8)), c("~~Futures-implied change", rep("",4), rep("Y",4)), c("~~Seasonal dummies", rep("",4), rep("Y",4)), 
           df = FALSE, digits = 2 , column.sep.width="0pt", omit.table.layout=c("nd"),   omit.stat =  c("F", "adj.rsq", "ser"), omit="Constant"  , star.cutoffs = c(0,0,0), 
           covariate.labels = c( "$ \\Delta \\log( \\frac{\\text{T-bill}_t}{\\text{Deposits}_t} )  $", "	$ \\Delta \\log( \\frac{\\text{T-bill}_{t-1}}{\\text{Deposits}_{t-1}} ) $", "$ \\Delta \\text{FFR}_t $", "$ \\Delta \\text{VIX}_t $ " ), 
           dep.var.caption = "\\textit{Dependent variable: $\\Delta$Treasury liquidity premium$_t$}", no.space=TRUE,
           order = paste0("^",vars.order,"$") )

par(mfrow=c(1,2))
MyPlot( TB_repoT$YearMon, TB_repoT$LiqPrem_Log, LegendNames = "log(liquidity premium)" )
MyPlot( TB_repoT$YearMon, TB_repoT$LiqPrem, LegendNames = "liquidity premium" )

par(mfrow=c(1,2))
MyPlot( TB_repoT$YearMon, TB_repoT$LiqPrem_LogDiff, LegendNames = "Diff log(liquidity premium)" )
MyPlot( TB_repoT$YearMon, TB_repoT$LiqPrem_Diff, LegendNames = "Diff liquidity premium" )


# ===================== Define TB_preGMM and GMM Functions ============================
# Critical: later analysis all rely on the TB_preGMM
TB_preGMM = TB[Is_Quarter_End==TRUE &  !is.na(LiqPrem) & !is.na(DebtToGDP) & !is.na(Deposits_GDP) & !is.na(VIX) & !is.na(FFR) ,  ]
TB_preGMM_KVJ = TB[ Is_Quarter_End==TRUE  & !is.na(LiqPrem) & !is.na(DebtToGDP) & !is.na(KVJ_fin_sector_debt_to_GDP) & !is.na(VIX) & !is.na(FFR) ,  ]
# View(TB[  , list(YearMon, LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ])

# ******** GMM Moment Function for the single equation  ******** 
g_single <- function(beta,X_matrix, weights=1, prediction=FALSE, use_debt_GDP_IV=FALSE, shut_off_VIX=FALSE, multiplicative_residual=FALSE, use_alternative_residual=FALSE, use_AR1=FALSE)
{
  X_matrix = as.matrix(X_matrix)
  rho = beta[1]
  beta_lambda = beta[2]
  # The structure of X_matrix = [  liquidity_premium, Deposit_spread, bond/GDP, deposit/GDP, consumption, GDP, VIX   ] 
  LiqPrem = X_matrix[,1]
  deposits_spread = X_matrix[,2]
  bonds_GDP = X_matrix[,3]
  deposits_GDP = X_matrix[,4]
  VIX = X_matrix[,5]
  if(use_debt_GDP_IV){
    total_debt_to_GDP = X_matrix[,6]  # Used as instrument
  }
  lambda = beta_lambda * VIX
  liq_prem_predicted = deposits_spread * (lambda*( bonds_GDP/deposits_GDP )^(rho-1) )
  residual = LiqPrem - liq_prem_predicted
  if(multiplicative_residual){
    residual = LiqPrem/liq_prem_predicted-1
  }
  if(use_alternative_residual){
    residual = LiqPrem/(lambda*( bonds_GDP/deposits_GDP )^(rho-1) ) - deposits_spread
  }
  if(use_AR1){
    u = LiqPrem/(lambda*( bonds_GDP/deposits_GDP )^(rho-1)) - deposits_spread
    # Assuming AR1 dynamics of the original residual u,  u_{t+1} = zeta*u_t + v_{t+1}, and v_{t+1} satisfies the orthogonality conditions.
    zeta = beta[3]
    residual = u - zeta*c( 0, u[1:(length(u)-1)] )
  }
  iv1 = VIX
  iv2 = deposits_spread
  if(use_debt_GDP_IV){
    iv3 = total_debt_to_GDP
    iv4 = total_debt_to_GDP^2
    mm = cbind( mm0=residual, mm_VIX=residual*iv1, mm_deposits_spread=residual*iv2, mm_totalTsy=residual*iv3, mm_totalTsy_square=residual*iv4 )
  }else{
    iv3 = bonds_GDP/deposits_GDP
    iv4 = NA
    mm = cbind( mm0=residual, mm_VIX=residual*iv1, mm_deposits_spread=residual*iv2, mm_bonds_to_deposits=residual*iv3 )
  }
  if(shut_off_VIX){
    column_list = 1:ncol(mm)
    mm = mm[ , column_list[-2] ]  # Remove the second moment with VIX
  }
  
  if(!prediction){
    return(mm*weights)
  }else{
    liq_prem_predicted_Only_FFR =   deposits_spread * mean( (lambda)*( bonds_GDP/deposits_GDP )^(rho-1) ) 
    liq_prem_predicted_Only_Supply = mean(deposits_spread) * (mean(lambda)*( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
    liq_prem_predicted_Only_VIX = mean(deposits_spread) * (lambda* mean(( bonds_GDP/deposits_GDP )^(rho-1)) )
    liq_prem_predicted_Only_Supply_and_VIX = mean(deposits_spread) * (lambda* ( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
    liq_prem_predicted_Only_Supply_and_FFR = deposits_spread * (mean(lambda)* ( bonds_GDP/mean(deposits_GDP) )^(rho-1) )
    liq_prem_predicted_Only_VIX_and_FFR = deposits_spread * ((lambda)* mean(( bonds_GDP/deposits_GDP )^(rho-1)) )
    LiqPrem_residual= LiqPrem-liq_prem_predicted
    deposits_spread_predicted = LiqPrem / (lambda*( bonds_GDP/deposits_GDP )^(rho-1) )
    return( cbind(data.table( residual, LiqPrem, liq_prem_predicted, deposits_spread, deposits_spread_predicted, LiqPrem_residual, iv1, iv2, iv3, iv4,  liq_prem_predicted_Only_FFR,liq_prem_predicted_Only_Supply,liq_prem_predicted_Only_VIX,liq_prem_predicted_Only_Supply_and_VIX,liq_prem_predicted_Only_Supply_and_FFR,liq_prem_predicted_Only_VIX_and_FFR), as.data.frame(mm)  ) )
  }
}


# ===================== TABLE 4: Core GMM Estimations =====================
beta_init = c( 0.3, 0.01 )
DT = TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX, DebtToGDP_raw ) ]
reg = gmm_customized(g_single, DT, beta_init, optfct = "nlminb", lower=c(0, 0), upper=c(30,1) )
result = g_single(reg$coefficients, DT, prediction = TRUE)
R2_vec = R2_fun( result$liq_prem_predicted, result$LiqPrem )

core_GMM_fun <- function(TBs, use_debt_GDP_IV=FALSE, shut_off_VIX=FALSE, beta_init=c(0.4,0.01), lags=4, use_AR1=FALSE, use_alternative_residual=FALSE, type="iterative"){
  N = length(TBs)
  R2_vec = rep(0, N)
  p_value_vec = rep(0, N)
  regs = list()
  for(iter in 1:N){
    DT = remove_missing(as.data.frame(TBs[[iter]]))
    reg = gmm_customized( function(beta,X) {g_single(beta, X, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX, use_AR1 = use_AR1, use_alternative_residual = use_alternative_residual)} , DT, beta_init, optfct = "nlminb", bw=lags, type=type )
    result = g_single(reg$coefficients, DT, prediction = TRUE, use_debt_GDP_IV = use_debt_GDP_IV, shut_off_VIX=shut_off_VIX, use_AR1 = use_AR1, use_alternative_residual = use_alternative_residual)
    R2_vec[iter] = R2_fun( result$liq_prem_predicted, result$LiqPrem )
    p_value_vec[iter] = round( as.numeric(summary(reg)$stest$test[2])  , digits=3  )
    regs[[iter]] = reg
  }
  return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}
# Regression Settings
TBs = list( TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ], 
            TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
            )
result = core_GMM_fun(TBs)
result_baseline = result
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)
result_with_IV_baseline = result_with_IV
  
result_core_list = c(result$regs, result_with_IV$regs)
reg_baseline = result_with_IV$regs[[1]]

#  ***** Main GMM Table *****
# Latex table
stargazer( c(result$regs, result_with_IV$regs), omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ), 
                   c( "Total Treasury IV?", "No", "No", "Yes", "Yes" ),
                   c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) , 
           dep.var.caption = "Measurement of B/D" )
# Print output
stargazer( c(result$regs, result_with_IV$regs) , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt",
          covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines =
            list( c("p-value of J-test",  c(result$p_value_vec , result_with_IV$p_value_vec) ),
                  c( "IV of $ \\frac{\\text{total debt}}{\\text{GDP}} $?", "No", "No", "Yes", "Yes" ),
                  c( "Variations explained", paste0(round( c(result$R2_vec, result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
          column.labels  = rep(c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),2) ,
          dep.var.caption = "Measurement of B/D", type="text" )





# ===================== TABLE 5: Explanatory Power of Different Components =====================
results = g_single( reg_baseline$coefficients, TB_preGMM[  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ], prediction=TRUE )
R2s = abs(c( R2_fun( results$liq_prem_predicted_Only_Supply, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_VIX, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_Supply_and_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted_Only_VIX_and_FFR, DT$LiqPrem  ),
R2_fun( results$liq_prem_predicted, DT$LiqPrem  )  ))
print( paste0("R2 = ", paste0(round(100*R2s,digits = 1),collapse="%, "), "%"  ), digits = 3 )


# ===================== Figure 4: Model Predictions versus Data =====================
result = g_single( reg$coefficients, DT[, c(1:5), with=FALSE ], prediction = TRUE )
# Panel A: Prediction on Liquidity premium
pdf( paste0(FigurePath,"Model_and_Data.pdf"), width = 5, height = 4 )
par(mar = c(2, 3.5, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, result[ , list(LiqPrem, liq_prem_predicted) ]  , ylab="Liquidity Premium (%)", LegendNames = c("Data", "Model Predictions") , LegendPosition = "topright" )
abline( v=as.Date(paste0(WWII_periods[c(1,length(WWII_periods))],"-01-01")), lty=2 )
text( as.numeric( as.Date( paste0(round(mean(WWII_periods)), "-01-01") ) ) , 2, labels = "WWII"   )
dev.off()

# Panel B: Prediction on FFR
pdf( paste0(FigurePath,"Model_and_Data2.pdf"), width = 5, height = 4 )
par(mar = c(2, 3.5, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, result[ , list(deposits_spread/deposits_spread_sensitivity_to_FFR, deposits_spread_predicted/deposits_spread_sensitivity_to_FFR) ]  , ylab="Federal Funds Rate (%)", LegendNames = c("Data", "Model Predictions") , LegendPosition = "topright" )
abline( v=as.Date(paste0(WWII_periods[c(1,length(WWII_periods))],"-01-01")), lty=2 )
text( as.numeric( as.Date( paste0(round(mean(WWII_periods)), "-01-01") ) ) , 20, labels = "WWII"   )
dev.off()



# ===================== TABLE 6: Subsample Analyses =====================
cutoff_date = min( TB_preGMM$YearMon[!is.na(TB_preGMM$Deposit_Spread_raw)] )
TBs = list( TB_preGMM[ !is.element(Year, WWII_periods)  , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[ !is.element(Year, WWII_periods)  , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[ Year>max(WWII_periods) , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[ Year>max(WWII_periods) , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[ Year>1980 , list( LiqPrem, Deposit_Spread, DebtToGDP, Deposits_GDP, VIX , DebtToGDP_raw) ],
            TB_preGMM[ Year>1980 , list( LiqPrem, Deposit_Spread, DebtToGDP, KVJ_fin_sector_debt_to_GDP, VIX , DebtToGDP_raw) ]
            )
result = core_GMM_fun(TBs)
result_with_IV = core_GMM_fun(TBs, use_debt_GDP_IV = TRUE)

stargazer( c(result_with_IV$regs) , omit.table.layout="n", star.cutoffs = c(0,0,0), df = FALSE, digits = 3 , column.sep.width="0pt", 
           covariate.labels = c( "$\\rho$", "$\\beta_\\lambda$" ), add.lines = 
             list( c("p-value of J-test",  c(result_with_IV$p_value_vec) ), 
                   c( "Variations explained", paste0(round( c(result_with_IV$R2_vec)*100, digits=1 ), "\\%") )  ),  dep.var.labels.include = FALSE,
           column.labels  = rep( c( "$ \\frac{\\text{Net Tsy}}{\\text{Deposits}} $", "$ \\frac{\\text{Net Tsy}}{\\text{KVJ Deposits}} $"),3) ,
           dep.var.caption = "Measurement of B/D" )



# ================ **SECTION 3** GMM Algorithm for the Nested Model =============================
g_inner <- function( beta, X_matrix, prediction=FALSE, use_debt_GDP_IV=FALSE ){
  X_matrix = as.matrix(X_matrix)
  eta = beta[1]
  beta_mu = beta[2]
  # The structure of X_matrix = [  liquidity_premium (P2_CP_Tbill), Tsy/GDP, VIX, KVJ_net_deposits(non_bank_debt_GDP),  LiqPrem_nonbank=P2_P1_spread ,   DebtToGDP_raw,  credit_spread, FFR ] 
  LiqPrem_Tsy = X_matrix[,1]
  Tsy_GDP = X_matrix[,2]
  VIX = X_matrix[,3]
  non_bank_debt_GDP = X_matrix[,4]
  LiqPrem_non_bank = X_matrix[,5]
  total_debt_to_GDP = X_matrix[,6]
  FFR = X_matrix[,7]
  
  LiqPrem_Tsy_predicted = beta_mu * VIX*(Tsy_GDP/non_bank_debt_GDP)^(eta-1)*LiqPrem_non_bank
  residual_inner = LiqPrem_Tsy - LiqPrem_Tsy_predicted
  
  iv1 = FFR
  iv2 = VIX
  if(use_debt_GDP_IV){
    iv3 = (total_debt_to_GDP)
    iv4 = (total_debt_to_GDP)^2
    mm = cbind( residual_inner, residual_inner*iv1, residual_inner*iv2, residual_inner*iv3, residual_inner*iv4) 
  }else{
    iv3 = Tsy_GDP/non_bank_debt_GDP
    mm = cbind( residual_inner, residual_inner*iv1, residual_inner*iv2, residual_inner*iv3) 
  }
  if(!prediction){
    return(mm)
  }else{
    return( data.table( LiqPrem_Tsy, LiqPrem_Tsy_predicted) )
  }
}

GMMs_inner_fun <- function(TBs, use_debt_GDP_IV=FALSE, rho_limit=2, remove_negative_liq_prem = FALSE, beta_init=c(0.5,0.1), type="twoStep"){
  N = length(TBs)
  R2_vec = rep(0, N)
  p_value_vec = rep(0, N)
  regs = list()
  for(iter in 1:N){
    DT = remove_missing(as.data.frame(TBs[[iter]]))
    if(remove_negative_liq_prem){
      negative_index = DT$LiqPrem<0 | DT$LiqPrem_nonbank<0
      print( paste0( sum(negative_index), " observations removed due to negative liquidity premium values!" ) )
      DT = DT[!negative_index,]
    }
    reg = gmm_customized( function(beta,X) {g_inner(beta, X, use_debt_GDP_IV = use_debt_GDP_IV)} , DT, beta_init, optfct = "nlminb", lower=c(0, 0), upper=c(rho_limit,1), type = type )
    result = g_inner(reg$coefficients, DT, prediction = TRUE, use_debt_GDP_IV = use_debt_GDP_IV)
    R2_vec[iter] = cor( result$LiqPrem_Tsy_predicted, result$LiqPrem_Tsy )^2
    p_value_vec[iter] = round( summary(reg)$stest$test[2], digits=3  )
    regs[[iter]] = reg
  }
  return( list( regs=regs, R2_vec=R2_vec, p_value_vec=p_value_vec )  )
}


sum(  TB_preGMM$P2CP3M -TB_preGMM$MMF_rate<0, na.rm=TRUE )

par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( AACP3M, P2CP3M ) ], LegendNames = c( "AA CP", "P2 CP" ), NA_elimination = FALSE )

pdf( paste0(FigurePath, "comparison_liq_prem.pdf"), width=8, height =4 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,2))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list(FFR - MMF_rate, Repo3M - MMF_rate, AACP3M - MMF_rate, P2CP3M - MMF_rate  ) ], LegendNames = c( "FFR-MMF rate", "Repo 3M-MMF rate", "AA CP-MMF rate", "P2 CP-MMF rate" ) )
abline(h=0)
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list(FFR - AACP3M, Repo3M - AACP3M, P2CP3M - AACP3M ) ], LegendNames = c( "FFR-AA CP rate", "Repo 3M-AA CP rate", "P2 CP-AA CP rate" ) )
abline(h=0)
dev.off()

pdf( paste0(FigurePath, "comparison_non_bank_volume.pdf"), width=5, height =4 )
par(mar = c(2, 2, 1, 1), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( MMF_to_GDP, KVJ_fin_sector_debt_to_GDP - Deposits_GDP ) ], LegendNames = c( "MMF/GDP", "(KVJ fin debt - deposits)/GDP" )  )
abline(h=0)
dev.off()

# Regression Settings
TB_reg = TB_preGMM[Year>=1976]

# Fed data: include P1 CP
par(mfrow=c(1,1))
MyPlot( TB_reg$YearMon, TB_reg[ , list( P2_CP_Tbill, P2_P1_spread,  P2_CP_Repo) ], LegendNames = c( "P2 CP - Tsy", "P2 CP - AA CP", "P2 CP - Repo" ) )

# Plot the relationship
par(mar = c(4, 4, 3, 3), mfrow=c(1,2))
plot( log(TB_reg$DebtToGDP/TB_reg$KVJ_net_deposits), log(TB_reg$P2_CP_Tbill/TB_reg$P2_P1_spread), xlab="log(Net Tsy/Net Non-bank Deposits)", ylab="log(P2CP-Tbill/P2-P1)" , mgp = c(2.2, 0.8, 0))
plot( log(TB_reg$DebtToGDP/TB_reg$KVJ_net_deposits), log(TB_reg$P2_CP_Tbill/TB_reg$P2_CP_Repo), xlab="log(Net Tsy/Net Non-bank Deposits)", ylab="log(P2CP-Tbill/P2-Repo)" , mgp = c(2.2, 0.8, 0) )

# Plot the ratios
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( P2_CP_Tbill/P2_P1_spread ) ], LegendNames = c( "P2 CP-Tbill/P2-P1 spread" ) )

pdf( paste0(FigurePath, "/CP_measures.pdf"), width=5, height = 4)
par(mar = c(3.5, 3.5, 1.0, 1.0), mfrow=c(1,1))
MyPlot( TB_preGMM$YearMon, TB_preGMM[ , list( P2_CP_Tbill,  P2_P1_spread, P2_CP_Repo ) ], LegendNames = c( "P2 CP-Tsy 3M spread", "P2 CP-P1 CP 3M spread", "P2 CP-Repo 3M spread" ) )
abline(h=0)
dev.off()


# ===================== Table 7: GMM Estimation of Inner Layer of the Nested System =============================
TB_reg = TB_preGMM[Year>=1970]
TBs = list( TB_reg[  , list( P2_CP_Tbill,    DebtToGDP , VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_P1_spread ,   DebtToGDP_raw,  credit_spread, FFR) ],
            TB_reg[  , list( P2_CP_Tbill,    DebtToGDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_CP_MMF,     DebtToGDP_raw,  credit_spread, FFR) ]
)

GMM_fun = GMMs_inner_fun
rho_upper = 1
result = GMM_fun(TBs, rho_limit = rho_upper, beta_init = c(0.5,0.1), type="twoStep")
result_with_IV = GMM_fun(TBs,use_debt_GDP_IV = TRUE, rho_limit = rho_upper,beta_init = c(0.5,0.1), type="iterative")

column_labels = rep(c("P2-P1 CP", "P2CP-MMF"),2)
lag_input = 4
stargazer( c(result$regs, result_with_IV$regs) ,  se=RobustSE(c(result$regs, result_with_IV$regs),lag_input = lag_input),
           add.lines = list( c("Variation Explained",  paste0(round(100*c( result$R2_vec, result_with_IV$R2_vec ),digits=1),"\\%") ),
                             c( "Total Treasury IV?", rep("No",length(result$regs)), rep("Yes",length(result_with_IV$regs)) ) ),
           dep.var.labels = paste0("GMM (NeweyWest with lag=", lag_input,")"),
           dep.var.caption=c(""), star.cutoffs = c(0,0,0), covariate.labels= c("$\\eta$", "$\\beta_\\mu$"), column.labels = column_labels  )


# # Use P2P1 spread as the non-bank liquidity premium measure
result_inner_list = c( result$regs, result_with_IV$regs )
result_inner = result_with_IV$regs[[1]]
eta = result_inner$coefficients[1]
beta_mu = result_inner$coefficients[2]
TB_preGMM[ , mu :=beta_mu*VIX/(1+beta_mu*VIX)  ]
TB_preGMM[ , bonds_GDP_inner :=  ( (1-mu)*KVJ_net_deposits^eta + mu*DebtToGDP^eta )^(1/eta)  ]
TB_preGMM[ , LiqPrem_bond_inner := ( (1-mu)^(-1/(eta-1))*P2_P1_spread^(eta/(eta-1)) + mu^(-1/(eta-1))*P2_CP_Tbill^(eta/(eta-1)) )^((eta-1)/eta)  ]

par(mar = c(3.5, 3.5, 1.0, 1.0), mfrow=c(1,2))
MyPlot(TB_preGMM$YearMon,  TB_preGMM[ , list( bonds_GDP_inner,DebtToGDP ) ], LegendNames = c("Constructed Inner Bond/GDP", "Net Tsy/GDP")  )
MyPlot(TB_preGMM$YearMon,  TB_preGMM[ , list( LiqPrem_bond_inner,P2_CP_Tbill ) ], LegendNames = c("Constructed Inner Bond LiqPrem", "Tsy Liq Prem")  )

MyPlot( TB$YearMon,  TB[ , list(Deposits_GDP +  DebtToGDP + KVJ_net_deposits  ) ] )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(Deposits_GDP)  ) ], LegendNames = c("FFR", "Deposits/GDP"), ylab="Standardized values"  )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(DebtToGDP )  ) ], LegendNames = c("FFR", "Net Tsy/GDP"), ylab="Standardized values" )
MyPlot( TB[ Year>1950]$YearMon,  TB[ Year>1950 , list(standardize(FFR), standardize(KVJ_net_deposits )  ) ], LegendNames = c("FFR", "Net Tsy/GDP"), ylab="Standardized values" )

# ===================== Table 8: GMM Estimation of the Whole Nested System =============================
g_nested <- function(beta,X_matrix, weights=1, prediction=FALSE, total_Tsy_IV=TRUE, use_VIX_moment=TRUE, alternative_residual=FALSE)
{
  # The columns of X_matrix are: (1) Tsy liquidity premium, (2) deposits spread, (3) bonds/GDP, (4) deposits/GDP, (5) VIX, (6) MMF/GDP,  (7) LiqPrem_MMF,  (8) Total Debt/GDP, (9) FFR
  X_matrix = as.data.frame(X_matrix)
  rho = beta[1]
  beta_lambda = beta[2]
  eta = beta[3]
  beta_mu = beta[4]
  # The structure of X_matrix = [  liquidity_premium, Deposit_spread, Tsy/GDP, deposit/GDP, VIX,  nonbank_deposits_GDP,  LiqPrem_nonbank, total_debt_to_GDP,  FFR] 
  LiqPrem_Tsy = X_matrix[,1]
  deposits_spread = X_matrix[,2]
  Tsy_GDP = X_matrix[,3]
  deposits_GDP = X_matrix[,4]
  VIX = X_matrix[,5]
  nonbank_deposits_GDP = X_matrix[,6]
  LiqPrem_nonbank = X_matrix[,7]
  total_debt_to_GDP = X_matrix[,8]
  FFR = X_matrix[,9]
  mu_t = beta_mu*VIX/(1+beta_mu*VIX)  # Important: remmeber how mu_t is approxed

  LiqPrem_Tsy_predicted = mu_t/(1-mu_t)*(Tsy_GDP/nonbank_deposits_GDP)^(eta-1)*LiqPrem_nonbank
  residual_inner = LiqPrem_Tsy - LiqPrem_Tsy_predicted
  if(alternative_residual){
    residual_inner = LiqPrem_Tsy/(mu_t/(1-mu_t)*(Tsy_GDP/nonbank_deposits_GDP)^(eta-1)) - LiqPrem_nonbank
  }
  
  LiqPrem_bonds = ( (1-mu_t)^(-1/(eta-1))*LiqPrem_nonbank^(eta/(eta-1)) + mu_t^(-1/(eta-1))*LiqPrem_Tsy^(eta/(eta-1)) )^((eta-1)/eta)
  bonds_GDP = ( (1-mu_t)*nonbank_deposits_GDP^eta + mu_t*Tsy_GDP^eta )^(1/eta)
  LiqPrem_bonds_predicted = deposits_spread*beta_lambda*VIX*( bonds_GDP/deposits_GDP )^(rho-1) 
  residual_outter = LiqPrem_bonds - LiqPrem_bonds_predicted
  if(alternative_residual){
    residual_outter = LiqPrem_bonds/(beta_lambda*VIX*( bonds_GDP/deposits_GDP )^(rho-1)) -  deposits_spread
  }

  iv_FFR = FFR
  iv4 = total_debt_to_GDP
  iv5 = (total_debt_to_GDP)^2
  iv_q_inner = Tsy_GDP/nonbank_deposits_GDP
  iv_q_outer = bonds_GDP/deposits_GDP
  if(total_Tsy_IV){
    if(use_VIX_moment){
      mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*VIX, residual_inner*iv_q_inner, 
                  residual_outter, residual_outter*iv_FFR, residual_outter*VIX, residual_outter*iv4, residual_outter*iv5 )
    }else{
      mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank,  residual_inner*iv_q_inner, 
                  residual_outter, residual_outter*iv_FFR, residual_outter*iv4, residual_outter*iv5 )
    }
  }else{
    if(use_VIX_moment){
      mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*VIX, residual_inner*iv_q_inner, 
                  residual_outter, residual_outter*iv_FFR, residual_outter*VIX, residual_outter*iv_q_outer )
    }else{
      mm = cbind( residual_inner, residual_inner*LiqPrem_nonbank, residual_inner*iv_q_inner, 
                  residual_outter, residual_outter*iv_FFR,  residual_outter*iv_q_outer )
    }
  }
  if(!prediction){
    return(mm*weights)
  }else{
    return( data.table( LiqPrem_bonds, LiqPrem_bonds_predicted, LiqPrem_Tsy, LiqPrem_Tsy_predicted) )
  }
}

coef_estimations = regs[[1]]$coefficients
R2_fun_nested <- function( coef_estimations, DT, total_Tsy_IV=TRUE ){
  result = g_nested(coef_estimations, DT, prediction=TRUE, total_Tsy_IV=total_Tsy_IV)
  R2_Tsy = R2_fun( result$LiqPrem_Tsy_predicted, result$LiqPrem_Tsy )
  R2_bond = R2_fun(  result$LiqPrem_bonds_predicted , result$LiqPrem_bonds )
  return( data.frame(R2_Tsy=R2_Tsy, R2_bond=R2_bond) )
}
R2s_fun_nested <- function( regs_list,  DT_list,  total_Tsy_IV_vec,  liq_prem_IV_vec ){
  R2_Tsy_vec = rep(0, length(regs_list))
  R2_bond_vec = rep(0, length(regs_list))
  for( iter in 1:length(regs_list) ){
    coef_estimations = regs_list[[iter]]$coefficients
    result = g_nested(coef_estimations, DT_list[[iter]], prediction=TRUE, total_Tsy_IV=total_Tsy_IV_vec[iter])
    # MyPlot( 1:nrow(DT_list[[iter]]), result[,list(LiqPrem_Tsy,LiqPrem_Tsy_predicted)] )
    # MyPlot( 1:nrow(DT_list[[iter]]), result[,list(winsorize_down_fun(result$LiqPrem_Tsy_predicted), winsorize_down_fun(result$LiqPrem_Tsy))] )
    R2_Tsy_vec[iter] =  summary( lm(winsorize_down_fun(result$LiqPrem_Tsy)~ winsorize_down_fun(result$LiqPrem_Tsy_predicted)) )$r.squared
    R2_bond_vec[iter] =  summary( lm(winsorize_down_fun(result$LiqPrem_bonds)~ winsorize_down_fun(result$LiqPrem_bonds_predicted)) )$r.squared
  }
  return( data.frame( R2_Tsy=R2_Tsy_vec, R2_bond=R2_bond_vec )   )
}


TB_for_reg = TB_preGMM

# note: use the projection of P2CP3M-Deposit rate on FFR as the deposit spread measure. 
TB_for_reg[ , CP_deposit_spread:= CP_deposit_spread_projection ]  

# Basic plots
par( mfrow=c(1,1) )
DT = TB_for_reg[Year>=1970 , list( P2_CP_Tbill,  CP_deposit_spread , DebtToGDP, Deposits_GDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_P1_spread, DebtToGDP_raw, FFR, YearMon) ]
DT = remove_missing(as.data.frame(DT))
DT1 = TB_for_reg[Year>=1970&Is_Quarter_End , list( P2_CP_Tbill, CP_deposit_spread , DebtToGDP, Deposits_GDP, VIX, KVJ_net_deposits,  LiqPrem_nonbank=P2_CP_MMF, DebtToGDP_raw, FFR, YearMon) ]
DT1 = remove_missing(as.data.frame(DT1))
columns = c(1:9)
# N = length(beta_init_list)

counter=0
beta_init_list = list()
for(iter in 1:length(result_core_list)){
  for(iter2 in 1:length(result_inner_list)){
    if( sum(result_core_list[[iter]]$coefficients>1) + sum(result_inner_list[[iter2]]$coefficients>1) == 0 ){
      counter = counter+1
      beta_init_list[[counter]] = c( as.numeric(result_core_list[[iter]]$coefficients), as.numeric(result_inner_list[[iter2]]$coefficients) )
    }
  }
}

regs=list()
regs[[1]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = FALSE) , DT[,columns], beta_init_list,  optfct = "nlminb", type="iterative" )
regs[[2]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = FALSE) , DT1[,columns], beta_init_list, optfct = "nlminb", type="iterative" )
regs[[3]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = TRUE), DT[,columns], beta_init_list, optfct = "nlminb", type="twoStep" )
regs[[4]] = gmm_customized_multi_init( function(beta, X) g_nested(beta,X,total_Tsy_IV = TRUE), DT1[,columns], beta_init_list, optfct = "nlminb", type="iterative" )
TB_R2 = R2s_fun_nested( regs,  list(DT[,columns],DT1[,columns],DT[,columns],DT1[,columns]), c(FALSE,FALSE,TRUE,TRUE) )
reg_nested_baseline = regs[[4]]
  
# Output to Latex
stargazer(regs, covariate.labels = c( "$\\rho$", "$\\beta_{\\lambda}$", "$\\eta$", "$\\beta_{\\mu}$" ), 
          star.cutoffs = c(0,0,0), dep.var.labels.include = FALSE, 
          dep.var.caption = "Measurement of non-bank liquidity premium", 
          add.lines = list( c("Total Treasury IVs?", "No","No","Yes","Yes"), 
                            c( "$R^2$ for Tsy Liq Prem", round(TB_R2$R2_Tsy,digits=2) ),
                            c( "$R^2$ for Q' Liq Prem", round(TB_R2$R2_bond,digits=2) ) 
          ),
          column.labels=rep(c("P2CP$-$P1CP", "P2CP$-$MMF"),2)
)
stargazer(regs, type="text", covariate.labels = c( "rho", "beta(lambda)", "eta", "beta(mu)" ), star.cutoffs = c(0,0,0), 
          column.labels=rep(c("P2CP$-$P1CP", "P2CP$-$MMF"),2), dep.var.labels.include = FALSE, 
          dep.var.caption = "Non-bank liquidity premium measures", 
          add.lines = list( c("Total Treasury IVs?", "No","No","Yes","Yes"), 
                            c( "R2 for Tsy Liq Prem", round(TB_R2$R2_Tsy,digits=2) ),
                            c( "R2 for Q' Liq Prem", round(TB_R2$R2_bond,digits=2) ) 
          )
)


# ===================== Figure 5: Plot the lambda and mu series ===================
beta_lambda_avg = 0
beta_mu_avg = 0
for( iter in 1:length(regs) ){
  beta_lambda_avg = beta_lambda_avg + regs[[iter]]$coefficients[2]/length(regs)
  beta_mu_avg = beta_mu_avg + regs[[iter]]$coefficients[4]/length(regs)
}
TB_for_reg[ , lambda := beta_lambda_avg*VIX/( 1+beta_lambda_avg*VIX ) ]
TB_for_reg[ , mu := beta_mu_avg*VIX/( 1+beta_mu_avg*VIX ) ]

pdf( paste0(FigurePath, "lambda_and_mu_raw.pdf"), width=8, height = 4 )
par(mar = c(3.5, 3.5, 1, 1), mfrow=c(1,2))
index = TB_for_reg$Year>=1970
MyPlot( TB_for_reg[index]$YearMon, TB_for_reg[index]$lambda )
abline(h=0.5)
MyPlot( TB_for_reg[index]$YearMon, TB_for_reg[index]$mu )
abline(h=0.5)
dev.off()


pdf( paste0(FigurePath,"mu.pdf"), width = 7, height = 4.5 )
par(mar = c(2, 3.5, 0.5, 0.5), mfrow=c(1,1))
lambdas_measured=TB_for_reg[index]$mu
lambdas_smooothed = rollmean(lambdas_measured, 40, fill = "extend")
MyPlot( TB_for_reg[index]$YearMon, data.frame(lambdas_measured,lambdas_smooothed ), ylim=c(0,1), ylab = "Liquidity Share of Treasuries",  LegendNames = c(expression(paste("Measured ", mu[t])),"Moving Average (10 years)") )
abline(h=0.5, lty=2)
dev.off()


# ===================== Table 9: Money Demand Regressions =============================

rho = reg_nested_baseline$coefficients[1]
beta_lambda = reg_nested_baseline$coefficients[2]
lambda =  beta_lambda*TB$VIX / (1+beta_lambda*TB$VIX)
eta = 1
beta_mu = reg_nested_baseline$coefficients[4]
mu = beta_mu*TB$VIX / (1+beta_mu*TB$VIX)

# Plot the slope coefficient on money demand. 
TB[ , logsp:=log(Deposit_Spread/(1+FFR/100)) ]
TB[ , loglsp:=log(LiqPrem/(1+FFR/100)) ]
TB[ , logY:=log(GDP_Real) ] # Log real GDP
TB[ , logcu := log(CURRCIR/GDP_nominal*GDP_Real) ] # currency in circulation
TB[ , Q1 := ((1-lambda)*(Deposits_GDP*GDP_Real)^rho + lambda*(DebtToGDP*GDP_Real)^rho)^(1/rho) ]
TB[ , NB :=  KVJ_fin_sector_debt_to_GDP - Deposits_GDP]
TB[ , B_prime := ( (1-mu)*(NB*GDP_Real)^eta + mu*(DebtToGDP*GDP_Real)^eta )^(1/eta)]
TB[ , Q2 := ( (1-lambda)*(Deposits_GDP*GDP_Real)^rho + lambda*(B_prime)^rho )^(1/rho)]
TB[ , logQ :=  log(Q1)-logY ] 
TB[ , logQ2 := log(Q2)-logY ] 
TB[ , logB := log(DebtToGDP) ]   
TB[ , logm := log(CURRCIR/GDP_nominal*GDP_Real + Checking_GDP*GDP_Real)-logY ]  # Including both currency and checking in the money demand
TB[ , logMK := log(KVJ_fin_sector_debt_to_GDP*GDP_Real)-logY ]
TB[ , logMD := log(Deposits_GDP*GDP_Real)-logY ]
TB[ , logMZ := log(MZM_to_GDP*GDP_Real)-logY ]
TB[ , logvix := log(VIX) ]
TB[ , adjsp := logsp - (rho-1)*log(KVJ_fin_sector_debt_to_GDP*GDP_Real) ]
TB[ , spread := LiqPrem/(1+FFR/100) ]
TB[ , QoverY := Q1/GDP_Real ]
TB[ , Q2overY := Q2/GDP_Real ]
TB_for_reg = TB[seq(1,nrow(TB),by=3)]
regs = list()
regs[[1]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year<1980 ]  )
regs[[2]] = lm( logm ~ logsp + logY , data = TB_for_reg[Year>=1980]  )
regs[[3]] = lm(logQ ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[4]] = lm(logQ ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
regs[[5]] = lm(logQ2 ~ logsp + logY + logvix, data= TB_for_reg[Year<1980] )
regs[[6]] = lm(logQ2 ~ logsp + logY + logvix, data= TB_for_reg[Year>=1980] )
stargazer(regs, se=RobustSE(regs), type="text", omit.stat =  c("F", "adj.rsq", "ser"), omit="Constant" )
stargazer(regs, se=RobustSE(regs), covariate.labels = c( "$\\log( \\frac{\\delta i_t}{1+i_t})$", "$\\log(\\text{real GDP}_t)$", "$\\log(\\text{VIX}_t)$" ), 
          star.cutoffs = c(0,0,0), dep.var.labels.include = FALSE, omit.stat =  c("F", "adj.rsq", "ser"), 
          column.labels = c(rep( c("pre 1980", "post 1980"), 3 ))  )


# =====================  [IMPORTANT] Comparison of Nagel 2016 Regressions =====================
# ================ TABLE 10: Replication of Nagel (2016) Table III =====================
# Nagel(2016) Replication
index = TB_full$Year<=2011
regs=list()
regs[[1]] = lm(  100*LiqPrem  ~  FFR + VIX ,  data=TB_full[index]  )
regs[[2]] = lm(  100*LiqPrem  ~  VIX + log(DebtToGDP_raw) ,  data=TB_full[index]  )
regs[[3]] = lm(  100*LiqPrem  ~  FFR + VIX + log(DebtToGDP_raw) ,  data=TB_full[index]  )
regs[[4]] = lm(  100*LiqPrem  ~  FFR + VIX + log(DebtToGDP) ,  data=TB_full[index]  ) 
regs[[5]] = lm(  100*LiqPrem  ~  FFR + VIX + log(DebtToDeposits) ,  data=TB_full[index]  ) 
stargazer( regs, omit = c("f"), omit.stat = c("F", "adj.rsq", "ser"), omit.table.layout="n",
           se=RobustSE(regs), df = FALSE, digits = 2 , column.sep.width="20pt", star.cutoffs = c(0,0,0), 
           covariate.labels = c( "FFR$_t$", "VIX$_t $", "log($ \\frac{\\text{Total Debt}_t}{\\text{GDP}_t} $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{GDP}_t} $)", "log($ \\frac{\\text{Net Tsy}_t}{\\text{Deposits}_t} $)"),
           dep.var.labels.include = FALSE)



# ===================== TABLE 11: Non-Linear Interaction Effects =====================
index =  !is.element(TB_full$Year,  WWII_periods) | is.element(TB_full$Year,  WWII_periods)   # All the periods. 
regs=list()
regs[[1]] = lm(  100*LiqPrem  ~  FFR + VIX  ,  data=TB_full[index]  )
regs[[2]] = lm(  100*LiqPrem  ~  VIX + log(DebtToGDP) ,  data=TB_full[index]  )
regs[[3]] = lm(  100*LiqPrem  ~  FFR + VIX + log(DebtToGDP) ,  data=TB_full[index]  )
regs[[4]] = lm(  100*LiqPrem  ~  FFR + VIX + FFR_and_Log_DebtToGDP ,  data=TB_full[index]  )
regs[[5]] = lm(  100*LiqPrem  ~  FFR + VIX + log(DebtToGDP) + FFR_and_Log_DebtToGDP ,  data=TB_full[index]  )
regs[[6]] = lm(  100*LiqPrem  ~  VIX_and_FFR + VIX_and_Log_DebtToGDP_and_FFR ,  data=TB_full[index]  )
stargazer( regs, omit = c("f"), omit.stat = c("F", "adj.rsq", "ser"), omit.table.layout="n",
           se=RobustSE(regs), df = FALSE, digits = 2 , column.sep.width="0pt", star.cutoffs = c(0,0,0), 
           covariate.labels = c( "FFR$_t$", "VIX$_t $", "log($ \\frac{\\text{Net Tsy}_t}{\\text{GDP}_t} $)" , "FFR$_t$*log($ \\frac{\\text{Net Tsy}_t}{\\text{GDP}_t} $)", 
                                 "VIX$_t $*FFR$_t $ ", "VIX$ _t $*FFR$ _t $*log($ \\frac{\\text{Net Tsy}_t}{\\text{GDP}_t} $)" ), dep.var.labels.include = FALSE )

